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ABSTRACT 

We present a new study of the Rayleigh-Taylor unstable regime of accretion onto 
rotating magnetized stars in a set of high grid resolution three-dimensional (3D) 
magnetohydrodynamic (MHD) simulations performed in low-viscosity discs. We find 
that the boundary between the stable and unstable regimes is determined almost en¬ 
tirely by the fastness parameter uj s = Q^/Q K {Ym), where Q* is the angular veloc¬ 
ity of the star and f^A:(r m ) is the angular velocity of the Keplerian disc at the disc- 
magnetosphere boundary r = r m . We found that accretion is unstable if uj s < 0.6. 
Accretion through instabilities is present in stars with different magnetospheric sizes. 
However, only in stars with relatively small magnetospheres, < 7, do the un¬ 

stable tongues produce chaotic hot spots on the stellar surface and irregular light- 
curves. At even smaller values of the fastness parameter, uj s < 0.45, multiple irregular 
tongues merge, forming one or two ordered unstable tongues that rotate with the an¬ 
gular frequency of the inner disc. This transition occurs in stars with even smaller 
magnetospheres, r m /R ;* < 4.2. Most of our simulations were performed at a small tilt 
of the dipole magnetosphere, 0 = 5°, and a small viscosity parameter a = 0.02. Test 
simulations at higher a values show that many more cases become unstable, and the 
light-curves become even more irregular. Test simulations at larger tilts of the dipole 0 
show that instability is present, however, accretion in two funnel streams dominates if 
0 > 15°. The results of these simulations can be applied to accreting magnetized stars 
with relatively small magnetospheres: Classical T Tauri stars, accreting millisecond 
X-ray pulsars, and cataclysmics variables. 
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1 INTRODUCTION 


Magnetospheric accretion occurs in different types of stars, 
including Classical T Tauri stars (CTTSs) (e.g.,[Bouvier et al. 


et al. 


2007), magnetized cataclysmic variables (CVs) (e.g., Warner 


|1995[ Hellier 2001), and accreting millisecond X-ray 


pulsars (AMXPs)(e.g., van der K lis 2006). The dynamically- 
important magnetic field stops the accretion disc at some dis¬ 
tance from the star, r m (called the magnetospheric radius), 
and subsequently, the magnetic field governs the flow of mat¬ 
ter onto the star (e.g.JPringle & Ree s||1972[ |Ghosh & Lamb| 
|1978[|Campbell|1992{|Lovelace et al.|1995|). 

The disc-magnetosphere boundary is prone to the Kelvin- 
Helmholtz and magnetic Rayleigh-Taylor (RT) instabilities 
(e.g., |Chandrasekhar||1961| |Arons & Lea||1976| ). It was sug- 
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gested in earlier studies that the instabilities at the disc- 
magnetosphere boundary may only lead to the mixing of 
plasma with the field in the external layers of the magneto¬ 
sphere |Arons 8c Lea 1976). However, our global 3D MHD 
simulations show that the disc matter can deeply penetrate 
the magnetosphere in tongues which produce chaotic hot 
spots on the stellar surface and irregular light-curves. Sim¬ 
ulations show that matter may accrete in either the stable 
or the unstable regime. In the stable regime, matter is lifted 
above the magnetosphere and accretes in two ordered fun¬ 
nel streams, forming two ordered hot spots on the surface of 
the star ( [Romanova et al.|2 003 , 2004; [Kuikarni & Romanovaj 
2005 ). In the unstable regime, matter penetrates through the 
magnetosphere in several equatorial “tongues” due to the 
magnetic Rayleigh-Taylor instability and forms several chaotic 
spots on the surface of the star ([Kuikarni & Romanova|2008{ 

Rom anova et al.|2008[|Bachetti et al.|2010[|Kurosawa & Ro-| 

manova 2013 ). The light-curves associated with the hot spots 
are expected to be periodic in the stable regime, and irregu- 
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Figure 1 . The boundary between stable and unstable regimes of accretion in the parameter space of r m and r cor for two different misalignment 
angles of the dipole: 6 = 5° and 6 = 20°. Here, the units of length are given in stellar radii for user convenience. Left panel: boundary line for 
0 = 5°. Stable, chaotic unstable, and ordered unstable accretion cases are represented by red squares, blue triangles, and green x’s, respectively 
Right panel: the boundary between stable and unstable regimes for © = 20°. The dash-dot-dot line shows the boundary for 0 = 5°. 


lar in the unstable regime, with a typical time-scale of a few 
peaks per rotational period of the inner disc. 


Simulations performed locally in the 3D simulation box 
also show that an ordered magnetic field is not an obstacle 
for the magnetic Rayleigh-Taylor instability (e.g., |Rastatter &| 
|Schindler||1999t |Stone & Gardiner||2007a|b). These s imula- 
tions and earlier simulations by Wang & Robertson] ([1984 , 
1985) show that the Rayleigh-Taylor instability leads to the 
formation of small-scale waves and filaments that merge to 
form much larger filaments, which then deeply penetrate the 
magnetically-dominated, low-density regions. 


The possible importance of the model for understand¬ 
ing CTTSs as well as other magnetized stars (such as AMXPs 
and IPs) pushed us to reconsider the unstable regime with 
new, more advanced simulations at a higher grid resolution. 
Test simulations showed that accretion through instabilities 
appeared in many more cases, including those which were 
stable in simulations with the coarser grid. For example, in 
the former simulations, accretion was usually stable in the 
cases of a small a—parameter of viscosity ( Shaku ra & Sun-| 
|yaev|1973] ) . However, in the new simulations, many of these 
cases become unstable. We concluded that the coarse grid was 
a factor that suppressed or weakened instability in many of 
our earlier simulations, so we reconsidered the problem us¬ 
ing a higher grid resolution. We also tested the convergence 
of the code at several grids with increasing resolutions. 


One of the important questions is, what determines the 
boundary between stable and unstable regimes of accretion? 
In prior work, we derived an empirical boundary that was 
based on the accretion rate (which we varied using the 
a—parameter of viscosity, e.g. Romanova et al. 2008 ). This 
boundary was useful, although it has not been completely un¬ 
derstood. This prior work led to the important conclusion that 
models with higher values of a are more unstable. One of the 
goals of our current research is to derive a boundary between 
stable and unstable regimes at a low viscosity in the disk. 


Another important phenomenon that has been systemati¬ 
cally observed in the new simulations is the frequent presence 
of what we call the ordered unstable regime, where the un¬ 
stable tongues merge to form one or two “ordered” tongues 
that rotate with the frequency of the inner disc. The accre¬ 


tion in ordered unstable regime has been observed in ear¬ 
lier simulations (performed at a lower grid resolution), but 
only in the cases of very small magnetospheres and high 
a—parameters (Romanova & Kulkarni 2009). In recent sim¬ 
ulations at a higher grid resolution, we observed that this 
regime is present at a wider range of parameter values, in¬ 
cluding larger-sized magnetospheres and lower viscosities in 
the disc. 

In this paper, we focus on two main issues: (1) searching 
a new boundary between the stable and unstable regimes; 
(2) investigating the ordered unstable regime and finding a 
boundary between the ordered and chaotic unstable regimes. 

In Sec. 2 we discuss the theoretical background of the 
problem. In Sec. 3 we describe the numerical model used in 
our simulations. In Sec. 4 we search for the boundary between 
stable and unstable regimes of accretion. In Sec. 5 we describe 
the ordered and chaotic regimes of unstable accretion, and 
derive the boundary between them. In Sec. 6 we analyze in¬ 
stability and investigate the dependence of instability on the 
corotation radius and the a—parameter of viscosity. In Sec. 7 
we investigate the grid convergence of the model. In Sec. 8 
we discuss the properties of accretion through instabilities in 
the cases of larger tilts of the dipole, 0. In Sec. 9 we compare 
relativistic and non-relativistic cases. In Sec. 10 we discuss 
the applications to different magnetized stars. We conclude 
in Sec. 11. 


2 THEORETICAL BACKGROUND 

2.1 Disk-magnetosphere boundary. 

The magnetic field of the star truncates the accretion disk at 
a radius r m where the magnetic stress in the magnetosphere 
matches the matter stress in the disk (e.g., [Pringle & Rees 
[T9721 |Lamb et al.|1973| ): 

p-\-pv 2 — B 2 /8it, or Pi = 8tt(p + pv 2 )/B 2 — 1 , (1) 

where p, p, v and B are the local density, gas pressure, ve¬ 
locity and magnetic field, respectively. Here, Pi is the gen¬ 
eralized plasma parameter, which includes both thermal and 
ram pressure (Ro manova et al.||2 002). It is analogous to the 
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Model 

At 

r COT 

a 

0 

P* 

Pinst 

Comments 

MO.3cl.805aO.O2 

0.3 

1.8 

0.02 

5 

2.4 

1.5 

intermediate (unstable chaotic/stable) 

,u0.3c2.5©5a0.02 

0.3 

2.5 

0.02 

5 

3.9 

1.9 

unstable ordered 

MO.3c501OaO.O2 

0.3 

5 

0.02 

10 

11.2 

2.2 

unstable ordered/chaotic 

MO.5cl.505aO.O2 

0.5 

1.5 

0.02 

5 

1.8 

0.9 

unstable chaotic 

MO.5c205aO.O2 

0.5 

2 

0.02 

5 

2.8 

2.4 

unstable chaotic 

MO.5c305aO.O2 

0.5 

3 

0.02 

5 

5.2 

2.3, 2.7 

unstable ordered/chaotic 

MO.5c505aO.O2 

0.5 

5 

0.02 

5 

11.2 

3 

unstable ordered 

MO.5c305aO.l 

0.5 

3 

0.1 

5 

5.2 

1.6 

unstable chaotic 

MO.5c505aO.l 

0.5 

5 

0.1 

5 

11.2 

2.6 

unstable ordered 

MO.5c302OaO.O2 

0.5 

3 

0.02 

20 

5.2 

1.9, 2.7 

unstable chaotic/ordered 

Atlcl.805aO.O2 

1 

1.8 

0.02 

5 

2.4 

0.85 

intermediate (unstable chaotic/stable) 

Atlc205aO.O2 

1 

2 

0.02 

5 

2.8 

0.4-1.6 

intermediate (unstable chaotic/stable) 

Atlc2.505aO.O2 

1 

2.5 

0.02 

5 

3.9 

3 

unstable chaotic 

Atlc305aO.O2 

1 

3 

0.02 

5 

5.2 

1.9, 3.8 

unstable ordered 

Atlc5©5a0.02 

1 

5 

0.02 

5 

11.2 

1.8, 7.4 

unstable ordered 

Atlcl.505aO.l 

1 

1.5 

0.1 

5 

1.8 

0.9 

unstable chaotic 

Atlc305aO.l 

1 

3 

0.1 

5 

5.2 

2.7 

unstable chaotic 

M2c2@5a0.02 

2 

2 

0.02 

5 

5.2 

2.8 

stable 

M2c3©5a0.02 

2 

3 

0.02 

5 

5.2 

2.8 

unstable chaotic 


Table 1. Sample (representative) models for different values of parameters g, r cor , a, and 0. The period of the star P* and the periods of 
instabilities Pi ns t obtained in Fourier spectral analysis are shown. 


standard plasma parameter p = 8irp/B 2 , but takes into ac¬ 
count the ram pressure of the matter flow in the disk. Axisym- 
metric and global 3D MHD simulations show that the condi¬ 
tion f3\ — 1 is valuable for finding the magnetospheric radius 
Tm (e.g., |Kulkarni & Romanova|2 013). This radius r m corre¬ 
sponds to the innermost edge of the disk, where the density 
drops sharply towards the low-density magnetosphere. This 
formula can be used to find r m in many situations, including 
those that are non-stationary or where the magnetic field is 
more complex than the dipole one. Sometimes, the condition 
P — 8ttp/B 2 = 1 is used to find the magnetospheric radius 
(e.g., |Bessolaz et aL| 2008). This condition, however, yields a 
somewhat larger radius at which matter flows from the disk 
to the funnel stream. 

In theoretical studies (e.g., Pringle & Rees|[l972l |Lamb| 
|et al.|[T973] ) the magnetospheric radius was estimated from 
the balance between the largest components of the stresses, 
assuming a dipole field of the star and Keplerian orbital speed 
in the disk; this gives 

r m = k[ni/(M 2 GM ir )] 1/7 , k~l, (2) 


Another important radius is the corotation radius, r cor , 
where, by definition, the angular velocity of the star matches 
the Keplerian angular velocity of the disc: = yjGM*/rl or , 

from which we obtain r cor = [GM*/Q* 2 ] 1//3 . Both the mag¬ 
netospheric radius r m and corotation radius r cor are useful in 
the analysis of magnetized stars. 

Another convenient parameter that is often used while 
investigating the physics at the disk-magnetosphere bound¬ 
ary is the fastness parameter, which is the ratio between the 
angular velocity of the star and the Keplerian angular velocity 
at the disk-magnetosphere boundary r = r m (e.g., |Ghosh &| 
|Lamb|1978[|Ghosh|2007| ): 


n* 


(r* m ) 


or, 



(3) 


One can see that this parameter is also a simple function of 
the ratio between two main radii, r m and r cor . The fastness 
parameter u s is expected to be significant in determining the 
processes at the disk-magnetosphere boundary (e.g., |Ghosh| 


2007). In our study we use either u s or the ratio r m /r c 


where /i* = P*P 3 is the magnetic moment of the star with a 
surface field P*, M is the disk accretion rate, and M* and P* 
are the mass and radius of the star, respectively. 

Axisymmetric simulations by, e.g., |Long et "ah] ( 2005 ) 
found that eqs.p^and 2|give similar values of r m if k « 0.5. 
|Kulkarni & Romanova (2013 ) performed a series of 3D MHD 


simulations of accretion in the stable regime and compared 
theoretically-derived values of r m (using Eq.[2]) with those ob¬ 
tained from simulations (using condition pi — 1). They found 
that, in the range of magnetospheric radii r m ~ (2 — 5)/?*, the 
radius r m obtained from the simulations can be described by 
Eq.[2]if 0.55 < k < 0.72. Here, we performed similar compar¬ 
isons between numerically-derived magnetospheric radii and 
Eq. fusing numerical runs in the stable and unstable regimes 
(see AppendixjA]). 


2.2 Rayleigh-Taylor instability at the 
disk-magnetosphere boundary 

The disk-magnetsophere boundary is prone to Rayleigh- 
Taylor and Kelvin-Helmholtz instabilities. Detailed simula¬ 
tions of the boundary in the shearing box (Rastatter & 
Schindler 1999) and also global simulations (e.g., Kulkarni & 
Romanova 2008) show that the RT instability strongly dom¬ 
inates. In a simplified case of two adjacent, rotating inviscid 
fluids with the magnetic field perpendicular to the boundary 
between them, the boundary is unstable if the effective grav¬ 
ity g eS — g + g c is negative (|Chandrasekhar| [l96l]). Here, 
g = —GM*/r 2 and g c = fl* r are gravitational and cen¬ 
trifugal acceleration, respectively. Taking into account that at 
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Figure 2. Left panel: 3D view of matter flow in a case where chaotic accretion in multiple tongues dominates, model ^lc2.505aO.O2, at time 
t = 19. One of the density levels is shown in color, selected magnetic field lines are shown in red. Middle panel: Same but in the face-on projection. 
Right panel: An equatorial slice of density distribution is shown in color. 
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Figure 3. An example of accretion in the chaotic unstable regime in the model /ilc2.5©5ct0.02 (times t = 19 — 21), where multiple tongues 
form. Top row shows consecutive 3D views of matter flow, where the color background represents one of the density levels and the lines are 
sample magnetic field lines. The axes show the directions for the angular momentum and magnetic momentum of the star. The middle panel 
shows consecutive xy slices. The color background shows density distribution and the lines show where the kinetic plasma parameter (3\ = 1. 
Bottom row shows the light-curve from rotating spots calculated at an inclination angle i = 45° (left panel), Fourier transform (middle panel) 
and wavelet transform (right panel) obtained from analysis of the light-curve. 


the boundary (r = r m ) g = -GM*/r 2 — -fl 2 K r, where 
LIk — \JGM*/r^ is the Keplerian velocity at the disc- 
magnetosphere boundary, we obtain 

9 eft — — l o 2 ). (4) 


One can see that in this simplified case the boundary is unsta¬ 
ble, if uj s < 1 . 

A more detailed theoretical analysis shows that in a real¬ 
istic situation there are a number of factors that can suppress 
instability. One of these factors is viscosity in the disc. This vis¬ 
cosity may be associated with the turbulence in the disk. Al¬ 
ternatively, it can be numerical viscosity associated with the fi¬ 
nite size of the grid in numerical models. Both types of viscos¬ 
ity can suppress the short wavelength perturbations ( |Chan-| 
drasekhar 1961). On the other hand, the azimuthal compo¬ 
nent of the field can also suppress the short-wavelength 
perturbation modes ( |Chandrasekhar|1 961). In addition to the 
above factors, the differential rotation of matter in the disk, 


dQ/dr, tends to suppress the RT unstable modes (e.g., Kaisig 
let al.|1 992). 


Another important factor that can increase or decrease 
the strength of instability is the compression factor if be = 

| In 1 , which is the gradient of surface density per unit of 
magnetic flux. The growth rate of the RT instability is larger 
if the gradient of T,/B z increases fast enough in the direction 


of gravity ([Kaisig et al. 

1992; Lubow & Spruit 1995). 

Spruit et al. (1995 

) have performed a general analysis of 


disc stability in the thin disc approximation, taking the veloc¬ 
ity shear into account The disc has a surface density E and is 
threaded by a magnetic field with the vertical component B z . 
Their analytical criterion for the development of instability is: 


7i3E = ( — Qeff) 


d , E 
— In —— 
dr B z 


> 2 



= 


(5) 


One can see that large, negative values of effective gravi¬ 
tational acceleration g e & and large values of the compression 
factor if be = | In -J^ | lead to stronger instability, while 
the radial shear, 7 ^ = 2 (r^) 2 , is a factor that suppresses 
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Figure 4. Left panel: Same as Fig. [ 2 ] but for a case where ordered accretion in two tongues dominates, model ^0.5c3©5a0.02, at time t = 14. 
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Figure 5. Same as Fig. [ 3 ] but for accretion in the ordered unstable regime, model /i0.5c3©5a0.02 (times t = 25.5 — 27.5), where one or two 
ordered unstable tongues form. 


instability by smearing out the perturbations. This criterion is 
valuable for the analysis of the disk-magnetosphere boundary 
However, it does not take into account the possible suppres¬ 
sion of perturbation modes by the azimuthal field, B or by 
viscosity Global 3D MHD simulations are required to take into 
account all of these factors, provided that the grid resolution 
is fine enough to resolve (and not suppress) the small-scale 
perturbation modes. 


3 NUMERICAL MODEL 


We perform global 3D MHD simulations of disc accretion onto 
a rotating magnetized star. The model we use is the same 
as in our earlier 3D MHD simulations of stable and unstable 
regimes (e.g., Romano va et al.||2008[ |Kulkarni & Romanova| 


2008), which has been described in our earlier papers (e.g., 
Koldoba et al.|2 002"| |Romanova et aH2 003| . Hence, we will 


only describe it briefly here. 


3.1 Initial and boundary conditions 

Initial conditions. A star has a dipole magnetic moment /i, 
the axis of which makes an angle © with the star’s rotational 
axis . The rotational axes of the star and the accretion disc 
are aligned. A star is surrounded by an accretion disc and a 


corona. The disc is cold and dense, while the corona is hot 
and rarefied, and at the reference point (the inner edge of 
the disc in the disc plane) the temperature and density are 
T c = 100 Td, and p c — 0.01 pd, where subscripts ‘d’ and ‘c’ de¬ 
note the disc and the corona. Initially, the disc and corona are 
in rotational hydrodynamic equilibrium (see e.g., Romanova 
|et al.||2002| for details). The disc is relatively thin, with the 
half-thickness to radius ratio h/r & 0.1. 

Boundary conditions. At both the inner and outer boundaries, 
most of the variables Aj are taken to have free boundary con¬ 
ditions, dAj/dr — 0. The free boundary conditions on the 
hydrodynamic variables at the stellar surface mean that ac¬ 
creting gas can cross the surface of the star without creating 
a disturbance in the flow. These conditions neglect the com¬ 
plex physics of interaction between the accreting gas and the 
star[j The magnetic field is frozen onto the surface of the 
star. That is, the normal component of the field, B n , is fixed. 
The other components of the magnetic field vary. At the outer 
boundary, matter is not permitted to flow into the region. The 
simulation region is usually large enough for the disc to have 


1 The impact of accreting matter with the stellar surface and forma¬ 
tion of the radiative shock wave near the surface has been investi¬ 
gated, e.g., by Koldoba et al. ( 2008| >; Orlando et al. ( |2Q1Q[ ). 
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Figure 6. The diagram shows transitions between stable and two types of unstable regimes discussed in the paper and calculated for low viscosity, 
a = 0.02, and a small tilt of the dipole magnetosphere, 0 = 5°. The diagram shows the dependence of boundaries on the fastness parameter lj 8 - 
The ratio rm/R-k shows the dimensionless size of the magnetosphere at which chaotic or ordered type of instability have been observed. 


enough mass to sustain accretion flow for the duration of the 
simulation run. 


3.2 MHD equations, numerical method and the grid 

The full set of 3D MHD equations taken in dimensionless 
form is solved numerically using a Godunov-type numerical 
scheme, written in a “cubed-sphere” coordinate system which 
rotates with the star (Kold oba et al.||2002] ). The numerical 
approach is similar to that described in |Powell| ( [1999| ), where 
the seven wave Roe-type approximate Riemann solver is used. 
The energy equation is written in the form of entropy balance, 
and the equation of state is that of an ideal gas. Compared 
with Kol doba et al.| ( |2002| ), viscosity terms are incorporated 
into the equations. Viscosity is modelled using the a-model 
( [Shakura & Sunyaev[ l973), and is incorporated only into the 
disc, so that it controls the accretion rate through the disc. In 
contrast with our earlier studies, we use a small a-parameter 
a = 0.02 in most of the simulation runs, and use a larger 
value a = 0.1 in the test runs. 

A “cubed sphere” grid consists of N r concentric spheres, 
where each sphere represents an inflated cube. Each of the 
six sides of the inflated cube has an N x N curvilinear grids 
which represent a projection of the Cartesian grid onto the 
sphere. The entire grid consists of 6N 2 x N r cells. The typical 
grid used in our simulations has N r — 140 cells in the radial 
direction, and N 2 = 61 2 ‘ angular’ cells in each block. This 
grid is twice as fine as the N 2 = 31 2 grid used in our earlier 
studies of instabilities (e.g., Kulkarni & Romanova 2008 ). To 
check the convergence of results at higher grid resolutions, 
the following grids were also calculated: 51 2 x 110, 71 2 x 
160, 81 2 x 185, and 101 2 x 220. Simulations show that at a 
small parameter of viscosity (a = 0.02) the coarsest grid, 31 2 , 
suppresses instability, while at all finer grids the results do not 
depend on the grid resolution (see details in Sec. 0 . 


3.3 Dimensionalization 

Equations are written using dimensionless variables. The di¬ 
mensionless value of every physical quantity Aj is defined as 
Aj = Aj/Ajo, where Ajo is the reference value for Aj. The 
simulations were performed in dimensionless variables A. 

First, we define the main reference values: the reference 
scale, Ro = R+/ 0.35, where R* is the radius of the star[^| the 
reference mass, M 0 = M*, where M* is the mass of the star; 
and the magnetic moment of the star, /z* = B*Rl, where B* 
is the magnetic field at the magnetic equator. Then, we de¬ 
termine the other reference values, which include: velocity 
vo — (GM 0 /Ro) 1/2 ; time-scale to = Ro/vo; period of rota¬ 
tion at r = Ro: Po = 2ttRq/vq ( used in our plots); refer¬ 
ence angular velocity £T 0 = v 0 /Ro ; and reference frequency 
fo = £To/27r. We determine the reference magnetic field B 0 
and magnetic moment po = B 0 Ro such that /i 0 = p*/jl, 
where p is the dimensionless parameter used to vary the size 
of the dimensionless magnetosphere. The reference field is 
B 0 = /x*/ (Rojd) , the reference density is po = Bq/vq, the ref¬ 
erence mass accretion rate is M 0 = povoRo = rI/(r 2 RoVo), 
and the reference surface density is S 0 = poRo . 

The results obtained in dimensionless variables A can be 
applied to different types of stars. We determine the dimen¬ 
sional mass, radius and magnetic field of a star and derive 
other reference values, as discussed above. Table [b] lists the 
reference values for three classes of stars: classical T Tauri 
stars, white dwarfs and neutron stars. To obtain the physi¬ 
cal dimensional values A, the dimensionless values A should 
be multiplied by the corresponding reference values A 0 as 
A = AA 0 . Subsequently, in the text we drop the tildes above 
the dimensionless variables and show dimensionless values 
everywhere unless otherwise specified. 


2 This reference scale has been used in our prior models (e.g., 
Koldoba et al. 2002; Romanova et al. 2002), and we now use it for 
consistency with our earlier works. 
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3.4 Main dimensionless parameters of the model. 

Taking into account the above dimensionalization, the num¬ 
ber of parameters in the problem is relatively small. The 
dimensionless mass and radius of the star, M* = 1 and 
R * 5 = 0.35, respectively, are fixed throughout all simulation 
runs. We fix the fiducial density values in the disk and corona 
at pd = 1 and p c — 0.01, respectively, and use the same disk 
structure across all simulation runs. The only parameter that 
is varied is the a—parameter of viscosity, which regulates the 
radial velocity in the disk, v r ~ a ( [Shakura & Sunyaevf l973 ). 
We use a small a—parameter, a = 0.02, in most simulation 
runs. 

There are two main dimensionless parameters that we 
vary. One of them is the dimensionless corotation radius r cor , 
which determines the dimensionless angular velocity of the 
star, O* = TcoJ 2 . This is one of the key parameters in deter¬ 
mining the boundary between stable and unstable regimes of 
accretion. 

The second important parameter is the dimensionless 
magnetic moment Jl, which determines the size of the dimen¬ 
sionless magnetosphere, r m = rm/Ro- We do not know the 
size of the magnetosphere in advance. Instead, we obtain it 
from the simulations. We can estimate the expected size of the 
magnetosphere using the theoretical formula Eq.[2]and substi¬ 
tuting the magnetic moment of the star, //*, with that obtained 
from our dimensionalization formalism: M 0 = pi / (Jl 2 Rqv 0 ), 
or pi = MqJI 2 RqVq and v 0 = (GM/Ro) 1 ^ 2 . We obtain: 


r m = rm/Ro = k(jf/M) 2/7 . 


( 6 ) 



Stable, r cor =2 Unstable, r cor =3 



Here, M is the dimensionless accretion rate, which we find 
from simulations and typically does not vary much across sim¬ 
ulations with similar a—parameters of viscosity. The main pa¬ 
rameter that determines the dimensionless size of the magne¬ 
tosphere, rm/R *, is parameter Jl. These two parameters (r cor 
and JT) determine the main physics at the disk-magnetosphere 
boundary [^] In the following sections, we drop the tilde’s 
above Jl and r cor for convenience. 

The results also depend on the tilt of the dipole magne¬ 
tosphere, 0. We investigated the cases of small (0 = 5°) and 
larger (0 > 20°) tilts in separate sets of runs. 


Figure 7. Top four panels: Example of two simulation runs with same 
parameters (p = 2, © = 5°, a = 0.02 at time t = 20) but different 
corotation radii (models p2c2@5a0.02 and p2c3G5a0.02). Left pan¬ 
els: xy-slice (top) and xz-slice (bottom) show the density distribution 
in case of r cor = 2, where accretion is stable. Right panels: accretion 
becomes unstable when r cor = 3. Red lines show sample magnetic 
field lines. Positions of the magnetospheric r m and corotation radii 
are shown with the dashed lines. Bottom four panels: Radial distribu¬ 
tion of of different terms of the Spruit criterion (Eq.[5]> in the vicinity 
of rm at time t = 6 (in the beginning of the unstable regime) for 
models shown above. 


4 BOUNDARY BETWEEN STABLE AND UNSTABLE 
REGIMES OF ACCRETION 

One of the main goals of present research was to find a set 
of parameters which would determine the boundary between 
stable and unstable regimes of accretion. 

4.1 Set of performed simulations 

To find the boundary between stable and unstable regimes of 
accretion, we performed a series of simulation runs, where we 
kept the same initial conditions for the density and pressure 
distribution in the disc and corona, but varied the the main 
dimensionless parameters of the model: magnetic moment p 
of the star (which determines the dimensionless size of the 

3 Note that the coefficient k equals to that in Eq. [ 2 ] 


magnetosphere) and the corotation radius r cor (which deter¬ 
mines the period of the star). We used a small parameter of 
viscosity, a = 0.02, in the majority of simulation runs to be 
sure that viscosity is not an essential factor in generating or 
suppressing instability. 

We ran two main sets of simulations for two different 
misalignment angles of the dipole moment: one at a rela¬ 
tively small angle, 0 = 5°, where instability is expected to be 
stronger, and the other at a larger angle, 0 = 20°, where fun¬ 
nel stream accretion is expected to be more favorable (e.g., 
|Kulkarni & Romanova|2 008). 

For each angle 0, we performed a series of simula¬ 
tions for stars with different corotation radii in the range 
of r co r = 1.2 — 5 and parameters p ranging from p = 0.1 
(small-sized magnetospheres) to p = 3 (large-sized magneto¬ 
spheres) . Table [l] shows sample models used for finding the 
boundary. Many of these models were also used for detailed 
analysis. The names of the models in the Table incorporate 
the parameters used in those models. 
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Figure 8. Left Panel: Dependence of different terms from Spruit’s criterion on corotation radius r cor . Simulations were performed for stars with 
fi = 1, a = 0.02 and different r cor . g e ft is multiplied by 10 for clarity Right Panel: Same but for dependence on a— parameter of viscosity 
Simulations were done for /x = 1, r co r = 1.5 and different a—parameters of viscosity 


To check the main hypothesis (which is the dependence 
of the boundary on the fastness parameter uj s and hence the 
ratio rm/rcor), we measured r m in each simulation run and 
marked whether accretion was stable or unstable at a given 
value of t C or (see Fig. [I). 

The accretion rate and r m generally vary in time. In some 
cases, accretion is either stable or unstable throughout the en¬ 
tire simulation run, while in other cases it is only marginally 
stable, and can transition from stable to unstable and back to 
stable again. In each case, we measured r m from the simula¬ 
tions using the position of the /5i = 1 line for a few moments 
in time, and marked whether the case was stable or unstable. 
The Pi = 1 line is not a smooth line: it reflects the inner parts 
of the disc as well as the matter-dominated unstable tongues. 
We determined r m using only the inner disc radius and ignor¬ 
ing the tongues. All of these points were used for finding the 
boundary between stable and unstable regimes of accretion. 


4.2 Results 

Fig. □ (left panel) shows the results of simulations in the 
r m — r cor parameter space for a small misalignment angle of 
the dipole, 0 = 5°. We can draw an approximate boundary 
line between stable (red squares) and unstable (blue triangles 
and green x’s) regimes of accretion. This line corresponds to 
the ratio r m /r cor ~ 0.71, or the fastness parameter u s « 0.6 
(see bold line in the figure). The plot shows that accretion is 
stable above this line, and is unstable otherwise. In this fig¬ 
ure, the values of r m and r cor are given in radii of the star for 
user convenience. As expected, instability occurs more easily 
when r m is smaller and r cor is larger, that is, when the mag¬ 
netosphere rotates more slowly than the inner disc and the 
fastness parameter co s is smaller. 

At a sufficiently small ratio of rm/r CO r, another transi¬ 
tion occurs and matter starts accreting in the ordered unsta¬ 
ble regime . The boundary between chaotic and ordered un¬ 
stable regimes approximately corresponds to the condition of 
rm jr cor ~ 0.59, or fastness parameter uj s ze 0.45. This new 
regime has only been observed in the cases of relatively small 
magnetospheres, r m /P* < 4.2. At larger magnetospheres, 
4.2 < Vm/R * < 7, accretion is chaotic (see left panel of 
Fig# 

We should note that at even larger magnetospheres, 
vm > 7 P*, the RT instability only develops in the external re¬ 
gions of the magnetosphere (e.g., |Romanova et al.|2014[|Ro-| 


Imanova & Qwocki[ 2015) and accretion in two funnel streams 
dominates. The parameter space with large r m should be in¬ 
vestigated separately 

We should also note that when r m /r co r > 1 (or, fastness 
uj s > 1), the magnetosphere rotates more rapidly than the 
inner disc, and the propeller regime is expected. We observed 
from the simulations that when r m « r COT (fastness parameter 
lj s = 1), the magnetosphere pushes the disc outward through 
the propeller mechanism. The parameter values under which 
the propeller regime is expected to dominate should be stud¬ 
ied separately. In the present study, we draw the uj s — 1 line 
as a reminder that the propeller regime is expected near or 
above this line. 


5 TWO TYPES OF UNSTABLE REGIME: CHAOTIC AND 
ORDERED 

Here, we describe the two types of unstable regime in greater 
detail, as observed in the simulations: (1) the chaotic unsta¬ 
ble regime, where matter accretes in several chaotic, transient 
tongues, and (2) the ordered unstable regime, where one or 
two tongues form and rotate persistently with the frequency 
of the inner disc. Below, we discuss these two types of unsta¬ 
ble regime for cases of a small tilt of the dipole, 0 = 5°. 

5.1 Chaotic unstable regime 

To demonstrate the chaotic unstable regime, we use the 
model /xlc2.505aO.O2, where the magnetosphere is relatively 
large {/i — 1) and the corotation radius is r cor = 2.5, which 
corresponds to the period P* « 3.9 in dimensionless units. 
We observed that matter accretes in several tongues, which 
form frequently and rapidly disappear. Fig. [2] (left two pan¬ 
els) shows a typical picture of matter flow through the chaotic 
unstable tongues, which are tall and narrow, and which pen¬ 
etrate the field lines by pushing them aside (right panel). 

Fig. HI shows the temporal evolution of tongues in sev¬ 
eral consecutive slices of density distribution, separated by 
the time interval of At = 0.5 Q Top two rows show that the 

4 Note that in our units At = 1 corresponds to the period of Keple- 
rian rotation at r = 1. However, the inner disc usually rotates with 
sub-Keplerian velocity, because the disc interacts with the slowly ro¬ 
tating magnetosphere 
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Figure 9. Top two rows: consecutive xy-slices (times t = 37 — 39) and analysis of variability for a larger misalignment angle, © = 10°, model 
^0.3c5©10a0.02. 3rd and 4th rows: same as top two rows, but for © = 20°, model /i0.5c5©20a0.02 (times t = 30.75 — 32.75). 


tongues are constantly modified, breaking or coalescing on a 
time-scale of At « 0.5 — 1. The 3D—panels show accretion 
through chaotic instabilities, with the occasional formation of 
funnel streams. 

The matter from unstable tongues accelerates due to 
gravity, falls onto the star and releases kinetic energy at the 
surface of the star, forming hot spots. We use a simple model 
for the radiation of the spots, suggesting that all kinetic en¬ 
ergy of the falling matter is converted into radiation, which 
is distributed isotropically (see details in Roman ova et al.| 
2004). We calculate the light-curve as seen by a remote ob¬ 
server, located at an angle of i — 45° with respect to the ro¬ 
tational axis of the star (see bottom left panel of Fig. [3]) . The 
light-curve shows irregular variations on a time-scale « 1.5 
times shorter than the period of the star, P* = 3.9. The bot¬ 
tom middle and right panels show the Fourier and wavelet 
spectra of this light-curve. The wavelet shows different angu¬ 
lar frequencies. One of them corresponds to the frequency of 
the star, /* = 0.26 (corresponding to the period of P* = 3.9). 
Other frequencies are higher, with /i ns t ~ 0.3,0.5. The Fourier 
spectrum shows that there are several peaks in the interval of 
Pinst = 2 — 3.2, which correspond to the wavelet frequencies. 
These periods also correspond to the time-scale of variability 
observed in the light-curve. We suggest that these periods re¬ 
flect the frequency of formation of the strongest tongues that 
reach the surface of the star. 

Both the wavelet and the Fourier analysis show the pres¬ 
ence of the stellar period, which can be associated with the 
fact that the entire set of unstable temporary spots is oriented 
around the magnetic pole. The magnetic pole is slightly tilted 
about the rotational axis (at © = 5° in this case) and the ro¬ 


tation of the star modulates the light from the set of chaotic 
spots. This is a probable reason for why the stellar period is 
also observed. Another possible reason is that some matter ac¬ 
cretes to the star in funnel streams above the magnetosphere 
and forms hot spots which tend to rotate with the star. In the 
present study, we chose an inclination angle of i — 45° be¬ 
cause this is where modulation by stellar rotation is expected 
to be the strongest. If this angle were smaller, then the stellar 
period would have had a smaller amplitude. 


5.2 Ordered unstable regime: one- and two-tongue 
accretion 

When the ratio r m /r co r becomes sufficiently small, unstable 
accretion becomes ordered: multiple tongues merge, forming 
one or two ordered tongues. In some cases, two tongues carry 
a comparable amount of matter flux, while in other cases one 
of the two tongues may carry more mass. In all cases, the 
bases of the tongues rotate with the frequency of the inner 
disc. 

Fig. 0 shows a snapshot of unstable accretion in two 
tongues, model /iO.5c305aO.O2. One can see that two ordered 
matter-dominated tongues push the magnetic field lines apart 
and penetrate into the deep layers of the magnetosphere, 
where they transition into regular funnel streams, but very 
close to the star. The top two panels of Fig. [5] show snapshots 
of accretion at five consecutive moments in time. The figure 
shows that initially there was only one tongue, and the sec¬ 
ond one formed later. One and two-tongue accretion modes 
are often observed in the same simulation run. 

The tongues deposit matter onto the stellar surface and 
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form one or two hot spots that rotate more rapidly than the 
star. The light-curve shows regular peaks that reflect the fre¬ 
quency of rotation of these spots (not the frequency of stellar 
rotation). The Fourier analysis shows a large peak associated 
with the rotating spots, with a period of Pi nst ~ 2 — 3. The pe¬ 
riod of the star, P * = 5.2, is also visible, but with a smaller am¬ 
plitude. The wavelet analysis shows frequencies in the range 
of /inst ~ 0.3 — 0.6, which correspond to the peaks observed 
in the Fourier spectrum. 

The main period observed in the light-curve is associated 
with the rotation of the ordered unstable tongues and approx¬ 
imately corresponds to the period of inner disc rotation. The 
ordered unstable regime may be an important mechanism in 
the formation of quasi-periodic oscillations during episodes of 
enhanced accretion, when the magnetosphere is compressed 
and the magnetospheric radius can be much smaller than the 
corotation radius. 


5.3 The boundary between chaotic and ordered 
unstable regimes 

To find the boundary between the chaotic and ordered un¬ 
stable regimes, we marked the cases in which one or two 
unstable tongues carry more mass than the other tongues 
and labelled these models as green x’s in Fig. [l] The ap¬ 
proximate position of this boundary corresponds to the ratio 
rm/r C or ~ 0.59 (the fastness o j s « 0.45 in the cases of rela¬ 
tively small magnetospheres, r m /P* <4.2. Ordered unstable 
accretion in one or two tongues dominates below this line. 
In the cases of larger-sized magnetospheres (r m /P* > 4.2), 
only the chaotic unstable regime has been observed. 

Simulations show that near the r m /r cor ~ 0.59 line ac¬ 
cretion is mainly chaotic, although at some point in a simu¬ 
lation run one or two tongues may start to carry more mat¬ 
ter than the other tongues. Accretion may also alternate be¬ 
tween ordered and chaotic within the same simulation run. 
At smaller values of r m /r CO r accretion becomes systematically 
more ordered, and the frequency associated with the ordered 
tongues becomes stronger. This means that at smaller values 
of the rm/r cor ratio the quality factor Q — //A/[^]of the 
quasi-periodic oscillations associated with the ordered unsta¬ 
ble tongues is higher. Such an increase in the quality factor 
is expected during accretion outbursts, when the disc moves 
inward towards the star. 

The diagram in Fig. [ 6 ] summarizes the boundaries be¬ 
tween the stable and unstable regimes, and between the 
chaotic and ordered unstable regimes. 


6 ANALYSIS OF INSTABILITY 

In this section, we analyze the dependence of the transition 
between stable and unstable regimes on different parameters, 
we are interested in understanding the dependence of accre¬ 
tion mode on the corotation radius r co r and the a—parameter 
of viscosity. 


5 The quality factor is the ratio between the frequency of oscillations 
and the width of the peak at the half-maximum. 


6.1 Comparison of stable and unstable cases 

To compare the simulation results with the theoretical cri¬ 
terion (Eq. we took two models with the same pa¬ 
rameters but different corotation radii: /x2c205c*O.O2 and 
/x2c305aO.O2, and compared their modes of accretion. In the 
first model, the corotation radius is r co r = 2 , and a star is 
in stable regime of accretion (see left panels of Fig. 0- In 
the second model, the corotation radius is r co r = 3, the star 
rotates more slowly, and the mode of accretion is unstable, 
(see right panels of Fig. [ 7 ]). We then calculated the different 
terms of the theoretical criterion (Eq. [5} in the vicinity of the 
magnetospheric radius r m , where the unstable perturbations 
can occur. To calculate these terms, we took the azimuthally- 
averaged values of relevant variables and plotted them as a 
function of radius in Fig. [ 7 ] Top panels show that in the stable 
regime (left panel, r cor = 2 ), 7 ^ > 7 ^ at the magneto¬ 
spheric radius (vertical dashed line), while in the unstable 
regime (right panel, r co r = 3), 7 > 7 ^. In both cases, 
results of the simulations are consistent with the analytical 
prediction (Eq. [ 5 } . Bottom panels show that in the unstable 
regime the effective gravity is g e & ~ —0.15, while in the stable 
regime it is slightly positive, g e ff ~ 0.06. 

The compression factor FTbe is larger in the unstable 
regime. Therefore, the term characterizing the instability, 
7 ^ s , is larger in the unstable regime due to both a larger 
effective gravity term — g e ff and a larger compression factor 
Pbe. 


6.2 Dependence of instability on r cor 

To better understand the physics of transition between the 
stable and unstable regimes and its dependence on the r cor 
parameter, we take the model /xlcl.505aO.O2 (which is at 
the boundary between stable and unstable regimes) and re¬ 
calculate it at different corotation radii r cor while keeping all 
the other parameter values fixed (/x = 1, 0 = 5°, a = 0 . 02 ). 
We then calculate the main terms from Eq. 1b e an d in, 
which determine whether the model is stable or unstable. 
We also separately calculate the compression factor, K B e = 

| In jU, and the effective gravity term, g eff . 

Fig. [ 8 ] (left panel) shows that effective gravity g e s de¬ 
creases with r CO v, changes sign from positive to negative, and 
then levels off at some negative value. It levels off because 
at large corotation radii the angular velocity is small and the 
centrifugal acceleration term in g e ff becomes negligibly small 
compared with the gravitational acceleration. One can also 
see that the compression factor K B e is large, but does vary 
systematically, and therefore the variation of the instability 
term 7^e is mainly determined by g e ff- We should note that, 
although both g e s and the shear term 7^ are relatively small 
and do not vary much with r cor , g e & becomes important due 
to its critical sign change during the transition. We conclude 
that the transition from the stable to the unstable regime is 
largely due to the variation in effective gravity. The transition 
occurs at r co r « 1 . 6 , when = 7 ^. 

We should note that the point of transition, r cor ~ 1.6, is 
close to the value of r co r ^1.5, where g e ff changes sign from 
slightly positive to increasingly more negative. 
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Figure 10. A comparison of simulation results for the model /rO.5cl.505aO.O2 (at time t = 18) using different grid resolutions. Simulations 
show that accretion is stable in the case of the coarsest grid, with 31 2 x 70 grids in one of the six sides of the cubed sphere (left panel ) and is 
unstable in the cases of the finer grids, 51 2 x 110, 61 2 x 140, 71 2 x 160, 81 2 x 185, and 101 2 x 220. The color background shows the density 
distribution in the equatorial plane. 


6.3 Dependence of instability on a 

In another set of experiments, we take the same model, 
/xlcl.505aO.O2, and vary parameter a in the range of a = 
0.02 — 0.1. We observed that this initially stable case becomes 
increasingly more unstable with increasing a. We are inter¬ 
ested in knowing why this transition happens and why the 
instability becomes stronger at higher values of a. 

Fig. [ 8 ] (right panel) shows results of our simulations. One 
can see that at larger a, the regime switches from stable to un¬ 
stable. The transition (where 7 ^ E = 7 ^) occurs at a « 0.037. 
At larger values of a, the disc-magnetosphere boundary be¬ 
comes more unstable due to both the decreasing effective 
gravity g eS and the increasing compression factor if be- The 
shear 7 ^ in the disc is relatively small and varies slowly com¬ 
pared with the 7 ^ e term, thus having a smaller influence on 
the transition. 

It is interesting to note that the compression factor iTss 
increases with a. This is probably due to the fact that the 
radial velocity of matter in the disc increases with a: v r = 
ac 2 s /v K (where c s and v K are the local sound speed and Ke- 
plerian velocity, respectively), and higher gradients of sur¬ 
face density (per unit of magnetic flux) are expected. On 
the other hand, the overall accretion rate increases with a: 
M ~ v r ~ a, which leads to a decrease in the magneto- 
spheric radius: r m ~ 1 /M 2/7 , and g e ff becomes more nega¬ 
tive. The transition between the stable and unstable regimes 
(where 7 ^ E = 7 ^) occurs at a « 0.037, which is near the 
point a = 0.028, where the effective gravity g e ff changes the 
sign. Here, as in the case of varying r co r , effective gravity plays 
an important role in the transition between stable and unsta¬ 
ble regimes. 

This set of simulations is similar to the simulations per¬ 
formed in our earlier studies, where we investigated the tran¬ 
sition between stable and unstable regimes by changing pa¬ 
rameter a to vary the accretion rate M, which we took to be 
the main factor in determining the mode of accretion (e.g., 
|Kulkarni & Romanova||2008| |Romanova et al.||2008| ). In this 
paper, our approaches were aimed at a deeper understanding 
of the physics of transition between the stable and unstable 
regimes, and its dependence on different parameters. 

In realistic discs, the angular momentum transport is 
probably provided by the magnetic turbulence, supported by 
the magneto-rotational instability ( [Balbus & Hawley| |l991). 
The effective a—parameter is determined by the ratio of mag¬ 
netic stress to matter pressure. The a— parameter may be rel¬ 


atively high in the inner parts of the disc, where the magnetic 
field of the turbulent cells is amplified by the rapid Keplerian 
rotation of inner disc matter (e.g., Hawley 2000; Armitage 
2002 ). In addition, part of the stellar magnetic flux penetrates 
into the disc and increases the magnetic stress in the inner 
disc, so that the a—parameter can be as high as a = 0.3 — 1 
( [Romanova et al.|2011||2012[ |Lii et al.|2014| ) . This is why we 
also performed test simulation runs at a relatively high pa¬ 
rameter of viscosity, a = 0.1. These simulations have shown 
that instability is stronger when a = 0.1 than in the cases of 
a = 0 . 02 , and that the boundary between stable and unstable 
regimes is expected to be higher, closer to the r m /r cor ~ 1 
line, though additional research is required to find this line. 


7 DEPENDENCE OF INSTABILITY ON GRID 
RESOLUTION 

To investigate the dependence of instability on grid resolu¬ 
tion, we chose a model with such parameters as make it close 
to the boundary between stable and unstable regimes of ac¬ 
cretion: /xO.5cl.505aO.O2. First, we ran simulations of this 
model using two grid resolutions, 31 2 x 70 and 61 2 x 140. We 
noticed that the 31 2 x 70 case is stable, while the 61 2 x 140 case 
is unstable. We then performed a set of simulations of this 
model at several other grid resolutions: 51 2 x 110, 71 2 x 160, 
81 2 x 185, and 101 2 x 220. We observed that instability is 
present at all finer grid resolutions, including 51 2 x 110 and 
finer (see Fig.fTo]). We compared matter fluxes at all grid res¬ 
olutions, and found that while matter flux values are sub¬ 
stantially different when comparing the 31 2 x 70 and the 
61 2 x 140 cases, there is smaller difference in the matter flux 
values of all the finer grid resolutions: 51 2 x 110, 61 2 x 140, 
71 2 x 160, 81 2 x 185, and 101 2 x 220 (see Fig. ]_1). These 
comparisons show that at a coarse grid resolution, numeri¬ 
cal viscosity suppresses the development of the RT instability, 
while at sufficiently fine grid resolutions, the results do not 
depend much on grid resolution. Numerical viscosity prob¬ 
ably plays the same role as regular viscosity in suppressing 
the smallest-scale perturbation modes, as discussed by |Chan-| 
|drasekhar| ( |1961| ). 


8 DEPENDENCE OF INSTABILITY ON 0 

The above simulations were performed for a small tilt of 
the dipole magnetosphere, 0 = 5°. It is important to know 
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whether accretion through RT instability is present in the 
cases of larger tilts. To answer this question, we performed 
a series of simulation runs for a larger tilt, © = 20°, with the 
same goal in mind: to find the boundary between stable and 
unstable regimes of accretion. Fig. [l] (right panel) shows the 
results of simulations. We observed that accretion becomes 
unstable in many instances, and the boundary between stable 
and unstable regimes corresponds to the ratio r m /r cor ~ 0.67 
(fastness parameter uj s ~ 0.54). This line is very close to the 
boundary line obtained for 0 = 5° (see dash-dot-dotted line 
right above the bold line in the right panel of Fig. [I]). This re¬ 
sult is in agreement with our hypothesis that g e ff is the main 
factor in determining the strength of the instability. 

To better understand the unstable accretion onto stars 
with different tilts of the dipole field, we performed test 
simulation runs at different tilts of the dipole: © = 

10°, 15°,20°,30°,40° and 60°. We focused on the ordered 
unstable regime and therefore used relatively small mag¬ 
netospheres, y — 0.3 - 0.5 and large values of r cor . One 
of the questions was whether this type of unstable accre¬ 
tion will take place in the cases of larger tilts of the dipole. 
First, we performed simulation runs for relatively small tilts, 
© = 10°, 15°, 20°, 30°. Fig.[9]shows xy —slices for two of these 
angles, © = 10° and 20°. The xy-slices of density distribution 
show ordered unstable accretion in one or two tongues in the 
case of the smallest tilt, © = 10° (top row of xy-slices). How¬ 
ever, at a larger tilt, © = 20°, accretion becomes more chaotic 
(bottom row of xy-slices). 

We also performed frequency analyses of the light-curves, 
which were calculated from the moving hot spots on the stel¬ 
lar surface, as seen at an inclination angle of i — 45°. The 
light-curve for the © = 10° case (see second row in Fig. 
shows that the main source of variability is the rotation of un¬ 
stable ordered spots. However, modulation by stellar rotation 
is also observed. The Fourier spectrum shows a peak with a 
period of Pi ns t ~ 2, associated with the instability, and a peak 
associated with the period of the star, P* « 11.2. The wavelet 
spectrum mainly shows the frequency /i ns t ~ 0.5, associated 
with the instability. 

In the case of © — 20° the amplitude of the oscillations 
associated with instability is smaller, and the light-curve is 
mainly determined by the rotation of the star, while the in¬ 
stabilities provide a high-frequency modulation of the main 
light-curve. The presence of instability is weak in the Fourier 
spectrum of the © = 20° case, but clearly visible in the 
wavelet spectrum (see bottom row of Fig. 12). 

In the cases of even larger tilts of the dipole, © = 30°, 
© = 40° and 60°, the light-curve is still modulated by the 
unstable component of accretion, but the amplitude of the 
modulation is smaller than in the case of 20°. 



Figure 11. Matter fluxes M onto the star as a function of time for 
model ^0.5cl.5©5a0.02, obtained at 6 different grid resolutions. 


in most of our prior simulations (Kulkarni & Romanova|2008, 


|2009t [Romanova et al.|2008]|Bachetti et al.|2010|). 

In this paper, we consider non-relativistic stars to be sure 
that the instability mechanism is present even in the cases 
of shallower, non-relativistic potential wells. However, in ap¬ 
plication to neutron stars it is important to take into account 
the relativistic effects and to compare the differences between 
instabilities in non-relativistic and relativistic cases. It would 
be too time-consuming to recalculate all the cases for rel¬ 
ativistic stars. Instead, we chose one representative model, 
/Hc3©5a0.02, and recalculated it taking into account the 
relativistic corrections. To model the relativistic effects we 
chose the Paczyriski-Wiita (PW) pseudo-relativistic potential 
( Paczyriski & Wiita 1980| ), <4>(r) = GM*/(r - r g ), where 
r g — 2 GM ~k / c z is the gravitational (Schwarzschild) radius. 
For a typical neutron star with mass M* — 1.4M 0 and radius 
P* = 10 km, the ratio r g /R * = 0.145. 

We found the results to be similar: instability of compara¬ 
ble strength has developed, with two ordered tongues form¬ 
ing and rotating with the frequency of the inner disc. FigEH 
shows xy-slices of these two models for the same moment in 
time. The matter flux onto the star (bottom panel) is about 
20% higher in the relativistic case. However, this does not af¬ 
fect the instability pattern qualitatively. We should note that 
the PW potential gives a steeper gravitational well than the 
fully relativistic approach, and therefore the relativistic effects 
will be even weaker in a more realistic relativistic case. 

We conclude that the main effect of enhanced gravity is 
enhanced accretion rate to the star, M. However, the differ¬ 
ence in accretion rate between the non-relativistic and the PW 
pseudo-relativistic cases is ~ 20%, which leads to a factor of 
only 1.2 1 / 7 « 1.03 difference in the magnetospheric radius 
r m . This difference is small, which is why relativistic effects 
do not enhance instability significantly. 


9 COMPARISON BETWEEN RELATIVISTIC AND 
NON-RELATIVISTIC CASES 

The results of current simulations are relevant for non- 
relativistic stars. However, in case of neutron stars, we 
should take into account the relativistic gravitational poten¬ 
tial, which provides a deeper gravitational well and may en¬ 
hance instability, in particular when the unstable tongues 
move closer to the star. A relativistic potential has been used 


10 APPLICATION TO DIFFERENT MAGNETIZED STARS 

The results of our simulations can be applied to different types 
of magnetized stars. 

10.1 Application to CTTSs 

Classical T Tauri stars stars show variability on different time- 
scales, ranging from hours to months ( [Herbst et al.||1994| ). 
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Figure 12. Top panels: density distribution of matter in the equatorial 
(xy) slice in the model with p = 1, r cor = 3 (model /ilc305aO.O2) at 
time t = 20 in non-relativistic (left) and relativistic (right) cases. Bot¬ 
tom panel: Matter fluxes to the star in non-relativistic and relativistic 
cases. 





10.2 Application to accreting white dwarfs 

A few types of accreting white dwarfs have strong magnetic 
fields. In the Intermediate Polars (hereafter IPs), the magne¬ 
tosphere disrupts the disc at r m > 10i?*, which is larger than 
the radii considered in this paper. Simulations of accretion 
to large magnetospheres (Romanova et aL]|2014| ) show that 
matter of the inner disc only partially penetrates the magne¬ 
tosphere, and most of the matter flows onto the star in two 
funnel streams (above the magnetosphere), forming two or¬ 
dered hot spots on the surface of the star. In these stars, the 
light-curves reflect the period of stellar rotation, and are ex¬ 
pected to be periodic. This process may explain the periodic 
light-curves of IPs. 

There is another class of accreting white dwarfs, called 
Dwarf Novae (DNs), that do not show period, but rather 
quasi-periodic oscillations are observed during accretion out¬ 
bursts (Warn er"et al.| |2004). The highest-frequency oscilla¬ 
tions (called Dwarf Novae Oscillations, or DNOs) have typ¬ 
ical time-scales of tens of seconds. The nature of these os¬ 
cillations has not yet been understood. It is possible that 
DNs are weakly magnetized white dwarfs whose magneto¬ 
spheres are strongly compressed during accretion outbursts, 
with v C or >> r m . This condition is favorable for ordered un¬ 
stable accretion, where one or two ordered tongues are the 
cause of quasi-periodic variability. Variation in accretion rate 
leads to variation in the position of the inner disc and fre¬ 
quency of the tongues, which may possibly explain the nature 
of DNOs. 


In most cases, the causes of variability have not been under¬ 
stood. 

Recent observations of the light-curves from young stars 
in the NCG 2264 cluster, obtained with the CoRoT telescope, 
show that many stars have irregular variability (Alencar et al. 
2010 ). One of the important time-scales corresponds to that 
expected in the unstable regime of accretion, that is, a few 
peaks per rotational period of the inner disc (|Stauffer et al.| 
|20l4l|Cody et al |2014^ . 

The ordered unstable regime may also be important dur¬ 
ing the periods of accretion outbursts, observed in EXor and 
FU Ori-type stars (e.g., |Audard et al.[ 2014 ). In these stars, the 
accretion rate is greatly enhanced, the magnetosphere is com¬ 
pressed, and the magnetospheric radius may be much smaller 
than the corotation radius. Therefore, conditions become fa¬ 
vorable for ordered unstable accretion in one or two ordered 
tongues. Recently, a very short period (of ~ 1 day) was dis¬ 
covered in the protostar VI647 Ori during its two recent ac¬ 
cretion outbursts ( [Hamaguchi et al.| 2012 ). Hamaguchi et al. 
(2012) suggested that the puzzling short-term X-ray variabil¬ 
ity of VI647 Ori was probably due to the rotational modula¬ 
tion of the hot spots on the stellar surface. However, as the 
authors noted, the ~ 1 day period corresponds to rotation at 
near-breakup speed of VI647 Ori. We suggest that it is im¬ 
probable that a star rotates so rapidly. Instead, the rapidly- 
rotating regions of enhanced plasma found in the empirical 
model of Hamaguc hi et al.| ( |2012| ) may be connected with the 
rotation of the ordered unstable tongues in the ordered unsta¬ 
ble regime of accretion, which is highly likely to occur during 
accretion outbursts. 


10.3 Application to accreting neutron stars 

The model can also be applied to accreting X-ray millisecond 
pulsars (AMXPs), where matter of the inner disk accretes onto 
a weakly magnetized (B ~ 10 8 — 10 9 G) neutron star during 
accretion outbursts (van der Klis|2006[|Ibragimov & Poutanenj 


2009] |Patruno & Watts|2012 ). Observations show pulsations 

associated with the rotation of the star, and also one or two 
peaks associated with the high-frequency quasi-periodic os¬ 
cillations. The frequencies of QPOs vary in time from values 
that are lower than the frequency of the star to values that are 
much higher than the frequency of the star. The origin of these 
QPOs is still not well understood. We suggest that, in AMXPs, 
one of the QPO frequencies may be associated with the rota¬ 
tion of unstable tongues in the unstable regime of accretion 
which onset when the accretion rate increases and the mag¬ 
netospheric radius approaches the value r m « 0.7r cor . The 
chaotic tongues rotate with the frequency of the inner disk, 
reflecting the QPO frequency (see also [Bachetti et al.] 2010 ). 
At even higher accretion rates, the disk moves even closer to 
the star and, at r m < 0.6r CO r, the ordered unstable regime 
dominates. The quality factor, Q, increases due the increas¬ 
ing dominance of the ordered unstable regime. This is consis¬ 
tent with the observations, which show that the frequencies of 
QPOs increase with X-ray flux (which is usually proportional 
to the accretion rate) (e.g., |Papitto et al.|2007| ), and the qual¬ 
ity factor increases with increasing QPO frequency (e.g., |van| 
|der Klis|2006l|Barret et al.|2007| ). 

We expect that the transition from the stable to the un¬ 
stable regime should lead to a decrease of the pulsed fraction 
in the observed radiation (e.g., |Kulkarni & Romanova| 2008). 
The pulsations are associated with stable (two funnel stream) 
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magnetospheric accretion, where two hot spots form and ro¬ 
tate with the frequency of the star. In contrast, in the unstable 
regime, most of the matter flows in irregular tongues, and 
the expected fraction of pulsed radiation is smaller. Interest¬ 
ingly, a drop in the pulse fraction has been recently observed 
in the AMXP SAX J1808.4-3658 ( |Bult & van der Klis|[2 015 ). 
The pulse fraction decreased dramatically when the QPO fre¬ 
quency (associated with the frequency of the inner disk) in¬ 
creased and became comparable to the stellar frequency. The 
drop in the pulse fraction may be connected with the tran¬ 
sition from stable to unstable accretion. However, these ob¬ 
servations show that the transition occurs when r m « r cor . 
This is different from the r m = 0.7r CO r boundary found in 
our simulations at a small parameter of viscosity (a = 0.02). 
Simulations show that at larger values of a (e.g., a — 0.1) 
the boundary corresponds to larger values of r m /r cor . Our 
model can explain the transition of the pulse fraction at high 
a—parameters of viscosity. 


11 CONCLUSIONS 

A new set of simulations has been performed with the goal 
of investigating the boundary between stable and unstable 
regimes of accretion and the properties of unstable accretion. 
Simulations were performed at twice as high a grid resolution 
as in our earlier studies ( [Kulkarni & Romanova|2008] , 2009; 
Romanova et al. 2008). A low viscosity parameter in the disc 

(a = 0.02) was used in most of the simulation runs. The main 
results of the new investigations are the following: 

1. We found that the boundary between stable and unsta¬ 
ble regimes of accretion depends almost entirely on the fast¬ 
ness parameter cj s (or the ratio between the magnetospheric 
radius and the corotation radius, r m /r cor ). Accretion is un¬ 
stable if (jj s < 0.6 (r m /r cor ) < 0.71), and is stable otherwise. 
The main simulations were performed at a small misalign¬ 
ment angle of the dipole, 0 = 5°, and low viscosity in the 
disc, a = 0.02. 

2. Below the lj s « 0.6 (r m /r cor ~ 0.71) line, accretion 
proceeds in the chaotic unstable regime through several unsta¬ 
ble tongues. The light-curve from the hot spots looks chaotic, 
and the spectral analysis shows several frequencies associated 
with the chaotic tongues. One of the frequencies often corre¬ 
sponds to the frequency of the inner disc, because the set of 
short-lived tongues rotates with the angular frequency of the 
inner disc. The frequency of the star is usually present in the 
frequency spectra if the simulation runs are sufficiently long. 

3. A new type of unstable regime of accretion, the or¬ 
dered unstable regime, was found in the cases of slowly rotat¬ 
ing stars, (jj s < 0.45 (r m /r cor < 0.59). In this regime, matter 
accretes in one or two ordered unstable tongues that rotate 
with the frequency of the inner disc. The ordered unstable 
regime has been observed in the cases of relatively small mag¬ 
netospheres, Vm/R * < 4.2. 

4. In the cases of larger misalignment angles of the 
dipole, accretion through instabilities is also present. The 
boundary between stable and unstable regimes was found for 
the case of 0 = 20°. This boundary, u s — 0.55, is very close to 
the one obtained in the case of a small tilt, 0 = 5°, because 
in both cases the instability is mainly determined by the ef¬ 
fective gravity. However, in cases of increasingly larger tilts, 
more matter flows above the magnetosphere in regular fun¬ 


nel streams and forms two hot spots near the magnetic poles, 
determining the regular sinusoidal light-curve that represents 
the period of the star. The matter accreting though instabili¬ 
ties provides the high-frequency modulation of the main light- 
curve. The amplitude of the oscillations decreases with 0 be¬ 
cause a strongly tilted dipole breaks the unstable tongues. 

5. At a higher viscosity in the disc, a = 0.1, chaotic in¬ 
stability becomes more irregular, and variability on different 
time-scales is observed in the light-curve. The frequency asso¬ 
ciated with the inner disc rotation is seen in both the Fourier 
and the wavelet spectra, while the frequency of the star has a 
much lower amplitude. 

6 . Analysis of the causes of instability in the borderline 
cases between stable and unstable regimes shows that: (a) 
Increasing the corotation radius r cor (while fixing all other 
parameters) leads to the larger negative values of effective 
gravity g e & causing the transition from the stable to the un¬ 
stable regime of accretion. We believe the sign change of # e ff 
to be the main factor in this transition, (b) Increasing the vis¬ 
cosity parameter a leads to a higher compression factor Abe 
which leads to stronger instability. 

7. A grid convergence has been observed at high grid 
resolutions. Simulations show that instabilities are only sup¬ 
pressed at the coarsest grid (31 * 1 2 3 4 angular cells in each block), 
while at all finer grids (51 2 and finer) accretion through in¬ 
stabilities is present, and is similar in all cases. 

8 . A comparison of relativistic and non-relativistic cases 
shows that instability is somewhat stronger in the relativistic 
cases. However, the relativistic potential is not the main factor 
in determining the mode of accretion. 

9. All of the above results were obtained for magneto¬ 
spheres with rm/R * < 7 and should not be generalized to 
stars with large magnetospheres. Separate studies show that 
in the cases of large magnetospheres, instabilities are only 
present in the external parts of the magnetosphere, while 
matter accretes to the star in two ordered funnel streams ( |Ro-| 
Imanova et al.|2014] ). 

10. The results can be applied to different types of accret¬ 
ing magnetized stars. The main findings are the following: 

(a) Accretion in ordered unstable tongues can be impor¬ 
tant during accretion outbursts, when the inner disc moves in¬ 
ward, compressing the magnetosphere, and the ratio r m /r cor 
is small. The frequency of the rotating tongues corresponds to 
the frequency of the inner disc and may be seen as a QPO fea¬ 
ture in the frequency spectra. The QPO frequency increases 
when the inner disc moves inward. The quality factor also 
increases during the inward motion of the disc. 

(b) When a star is observed for a few periods of stellar 
rotation, as in CTTSs, the frequency of the oscillations associ¬ 
ated with the ordered unstable tongues may be mistaken for 
the frequency of the star. 

(c) A star may alternate between stable and unstable 
regimes of accretion, thus showing an intermittency in its pul¬ 
sations. 
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APPENDIX A: MAGNETOSPHERIC RADIUS 

It is important to know the relationship between the magne- 
tospheric radii obtained from the balance of stresses, /3i = 1 
(which have been obtained from our 3D MHD numerical sim¬ 
ulations and used in our paper), and the magnetospheric radii 
derived from the theoretical formula, [2] which is widely used 
for estimating r m (in the cases where the values of /x*, M and 
M* are known). In our simulations, we can compare these 
two radii by (1) taking one of them from the condition /A = 1, 
and (2) taking the second one from Eq. [2] using parameters 
/x*, M* and M, whose values are known from our simula¬ 
tions. Comparisons of these radii will help us derive the un¬ 
known coefficient k, as well as check the overall dependencies 
in Eq. [2] 

For these comparisons, we first convert the theoretically- 
derived magnetospheric radius (Eq. [2]) to the dimensionless 
form (see procedure of dimensionalization in Sec. |3.4| and 
resulting Eq.[6])0 We re-write Eq.[6]in a more general form: 

Tm k(n 2 /M) n , (Al) 

and search for coefficients k and powers n. All variables are 
dimensionless, but we removed the tilde’s for convenience. 
Note that coefficient k is exactly the same as in Eq. [2] and the 
power n can be different from n — 2/7 used in the theoretical 
formula. 

Next, we take our models with different values of param¬ 
eters /i and t cor ? and find r m from the condition /3i = 1. We 
also measure the dimensionless matter flux M and calculate 
the value /i 2 /M. We then plot r m versus /i 2 /M. /i is a param¬ 
eter of the model, and is fixed in each model. However, the 
accretion rate M varies in time. In some cases, the accretion 
rate levels off at a constant value, in spite of the regime be¬ 
ing strongly unstable, as shown in the left panel of Fig. |A1| 
Alternatively, it may vary strongly (by a factor of 2-3), as in 
the cases of the chaotic unstable regime. This is why the main 
dispersion of points is expected to arise from the variation of 
M. We often took several different points in time within the 
same simulation run to see how the variation in accretion rate 
would affect the resulting radius 0 

To measure the magnetospheric radius, we used the con¬ 
dition (3i = 1. However, in the unstable regime, the magneto¬ 
spheric radius often has a complex shape, and the /3i = 1 line 
is not truly circular. To obtain r m we approximate the inner 
boundary using the procedure shown in Fig. |A1[ where the 
circle approximately reflects the inner boundary. We should 
note that, in any particular simulation run, the value of r m 
does not vary much, which may be the result of a weak (ex¬ 
pected) dependence on the accretion rate: r m ~ M~ 2 ^ . For 


example, variation of M by a factor of 3 leads to the variation 
of r m by a factor of 1.3. 

Figure [A2| shows the results of simulation runs in the sta¬ 
ble and unstable regimes (left and right panels, respectively). 
One can see that in both cases the magnetospheric radius 
increases with parameter /j 2 /M, and the best-fit curve is a 
power-law (see dashed lines), with the power of n = 0.22 in 
the unstable regime and n — 0.19 in the stable regime. The 
plot also helped derive coefficient k, which is k « 0.78 in the 
stable regime and k « 0.89 in the unstable regime. Therefore, 
we expect the formula for the magnetospheric radius to be in 
the following form: 


0.78 [/4/(M 2 GM*)] c 


(A2) 


0.891/4/(M 2 GM*)] c 


(A3) 


In both formulae, the power n is smaller than the power 
1/7 « 0.146 in theoretical formula [2] This issue has been 
analyzed in detail by Kulkarni & Romanova (20130 who 
concluded that the power n is smaller than in the theoreti¬ 
cal formula due to the compression of the magnetosphere by 
the disc and the non-dipole shape of the external regions of 
the magnetosphere. This may also be the reason for why in 
the stable regime the power n is smaller than in the unstable 
regime: in the unstable regime, the compression is expected 
to be smaller due to the penetration of matter in unstable 
tongues. 

We can estimate the difference between using the the¬ 
oretical formula [2 and formulae [A2] or |A3| For compar¬ 
isons, we take Eq. Al] in the form for stable regime, r m = 
0.78(/i 2 /M) 019 and in the form corresponding to Eq. [2] 
rm = k(ii 2 /M) 2/7 and find k = 0.78 (/tt 2 /M)~ 0 - 095 . Then 
we note, that in our model, varies in the range of 

2 — 6 which corresponds to 0.2 < / 1 2 /M < 80. For this 
range, we obtain 0.50 < k < 0.90. Similarly, for the un¬ 
stable regime we find k = 0.89(/i 2 /M) -0 026 and condition 
0-80 < k < 0.93. One can see that in the case of small mag¬ 
netospheres 2 < rm/R * < 6, theoretical formula [2] can be 
used to describe the stable regime, if 0.50 < k < 0.90, and 
the unstable regime, if 0.80 < k < 0.93. 


APPENDIX B: REFERENCE VALUES IN THE TABLE 

Table [b] shows sample values of physical parameters for dif¬ 
ferent types of stars. See Sec.[3]for a detailed description. 
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Figure Al. Left panel: Temporal variation of the dimensionless matter flux M in one of models in the unstable regime, ^0.5c2.5©5a0.02. Middle 
and right panels show the equatorial slices of the density distribution at time t — 15 and t = 21. The thin purple line corresponds to condition 
/3i = 1 (see Eq. [lj . Thick solid line shows the approximate position of the magnetopsheric radius. 





Figure A2. Left panel: The magnetospheric radii taken from simulations in the stable regime, r m (measured in stellar radius R *) versus the value 
p 2 /M. Right panel: same, but for the unstable regime. Dashed lines represent the best power fit. 
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Table Bl. Sample values of physical parameters for different types of 
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